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Abstract: This paper presents numerical methods for computing regions of finite-time 
invariance (funnels) around solutions of polynomial differential equations. First, we present a 
method which exactly certifies sufficient conditions for invariance despite relying on approximate 
trajectories from numerical integration. Our second method relaxes the constraints of the first by 
sampling in time. In applications, this can recover almost identical funnels but is much faster to 
compute. In both cases, funnels are verified using Sum-of-Squares programming to search over a 
family of time- varying polynomial Lyapunov functions. Initial candidate Lyapunov functions are 
constructed using the linearization about the trajectory, and associated time-varying Lyapunov 
and Riccati differential equations. The methods are compared on stabilized trajectories of a 
six-state model of a satellite. 

Keywords: Nonlinear systems, Lyapunov methods, Stability domains 
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1. INTRODUCTION 

In this paper we propose algorithms for computing a 
finite-time "funnel" of a dynamical system: a set of initial 
conditions whose solutions are guaranteed to enter a 
certain goal region at a particular time. This work is 
motivated by new algorithms for nonlinear control design 
wherein the controller is constructed from a tree of finite- 
time trajectories, each locally stabilized and all leading to 
a certain goal point (see Tedrake et al. (2010)). Sparseness 
of this tree is advantageous, and is directly related to the 
size of the funnel that can be verified for each trajectory 
in the tree. 

The basic method is to search over a class of time- varying 
Lyapunov functions about the trajectory, and verify a 
positive- invariance condition via Sum-of-Squares program- 
ming. A preliminary algorithm for this verification was 
suggested in Tedrake et al. (2010). This paper extends that 
work by replacing a greedy set of nonconvex optimizations 
by an alternation of convex optimizations (similar to that 
proposed for equilibria in Jarvis-Wloszek et al. (2005)). It 
is shown that the verified regions can be made exact even 
if the trajectory is an approximate (numerical) solution. 
Furthermore, an alternative time-sampled verification is 
suggested which recovers almost identical funnels and is 
substantially faster to compute. 

Recently, a great deal of research has been dedicated to 
region-of-attraction analysis on polynomial vector fields 
(Papachristodoulou and Prajna (2002); Parrilo (2003); 
Topcu et al. (2008); Tan and Packard (2008)) and more 
general non-polynomial vector fields (Papachristodoulou 
and Prajna (2005); Chesi (2009)) through Sum-of-Squares 
programming. There is comparatively little written in 
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the Sum-of-Squares literature addressing time- varying sys- 
tems. In Julius and Pappas (2009) finite time-invariance 
around trajectories is explored to provide outer approxi- 
mations of the set reached from some initial conditions. 
This is accomplished by using regionally valid Lyapunov 
certificates to construct barrier functions bounding expo- 
ncntia, certifying solutions do not enter keep-out sets. By 
contrast, the algorithm in this paper constructs inner- 
approximations of the solutions which can reach a goal 
region through the construction of time- varying Lyapunov 
functions. 

1.1 Notation 

We denote the set of n-by-n positive definite matrices 
by §™. For P e $1 and x e R" we use \\xf P as 
short hand for x'Px. We denote the set of polynomials 
in x € R™ with real coefficients by R[x]. The subset 
of these polynomials which are Sum-of-Squares (SOS) is 
denoted Tt[x], that is: p(x) e T,[x] if and only if there 

exists {<?» C R[x] such that p(x) = J2i=i 9i( x ) 2 - 
We similarly denote polynomials and SOS polynomials in 
multiple vector valued variables (x l7 . . . , x/.) € R™ 1 x . . . x 
R nk by R[xi, . . . , x k ] and Y,[x\, . . . ,Xk]- Finally, for a set 
A C R™, int(A) refers to its interior. 

2. PRELIMINARIES 
We are given a time-varying nonlinear dynamical system: 

|*(t) = /(*,*(*)), (1) 

where / : [t^, t{] x R" i->- R™ is piecewise polynomial in t 
and polynomial in x. Further we are given an bounded 
"goal region", Q C R™, with non-empty interior. In 
applications, we will generally be investigating systems (1) 



arising as closed loop systems stabilizing some trajectory. 
For convenience we make the following definitions. 

Definition 1. Given the dynamics (1), a set T C [io,if] x 
K™ is a funnel if for each (r, x T ) in T , the solution to (1) 
with x(t) — x T exists on [r, if] and for each t € [r, if] we 
have (t, x(t)) e J 7 . 

Definition 2. Given the dynamics (1) and the goal region 
Q, a set T C [to,t{] x K n is a funnel into Q if it is a funnel 
and for any x £ M" wc have (tf, x) € T implying x e Q. 

We see that a funnel into Q is an inner-approximation 
of the set of solutions which flow through Q at the time 
t { . In this work we are interested in finding the largest 
possible funnel into Q , measured, for example, by volume 
as a subset in [to, tf] x W 1 . Our approach uses time- varying 
Lyapunov functions, exploiting the following Lemma. 

Lemma 3. Let V(t,x) be a function V : [t ,tf] x R™ >->• 
[0, oo) piecewise continuously diffcrentiable 1 with respect 
to t and continuously differentiable with respect to x. For 
each t € [to,t{] we define: 

fi t := {x | K(t,x) < 1}, 

and 

dfi* := {x | V(t,a;) = 1}. 
If for each t e [to, if] and a; G 9fi t we have: 

— V(t,x)f(t,x) + -V(t,x)<0, (2) 

then the set: 

F={(t,x)\x€fl t } (3) 
is a funnel. If additionally fit f C then J 7 is a funnel into 

To leverage the SOS relaxation, our goal region Q will be 
the closed interior of an ellipse: 

Q = I Nik < i}, 

where Pq E M. nxn is a symmetric positive definite matrix. 
More general goal regions can be defined as sub-level sets 
of polynomials. 

2.1 A Class of Lyapunov Functions 

We begin by describing a class of candidate Lyapunov 
functions based on solutions of (1) and related Lyapunov 
differential equations. To leverage SOS we approximate 
these solutions by piecewise polynomials in time. This 
approximation does not render our certificates inexact. 
Further, we describe how sufficiently fine approximation 
of these solutions guarantees the existence of a certificate. 

We assume we have access to a nominal solution of (1), 
x ■ [t, if] i ^ R n with r e [to, if), such that Xf := x (tf) e 
int (<"/). We solve for a symmetric positive definite matrix 
Pf G S™ describing the largest volume ellipse centered on 

Xf, 

£ { = {x | \\x-xf\\% < 1}, 

contained in the goal region (i.e. £f C Q). When Q is an 
ellipse, this containment problem is an SDP. Given a more 
general polynomial goal region, one can relax the problem 

1 We will generally write conditions on JrV. At points of disconti- 
nuity, wc require these conditions to hold for both the left and right 
derivative of V with respect to t. 



to finding the largest sphere centered on Xf contained in Q 
to a SOS program. 

Our Lyapunov functions are parameterized by a time- 
varying, symmetric positive definite matrix P* : [r, if] 
§":' 

F*(t, a; ) = || a; - a :o(t)||^ (t) . (4) 

We require P*(tf) > Pf, so that fi tf , the one sub-level set 
of V*(tf, x), is a subset of Ef and thus contained in the goal 
region. 

Our optimization approach requires an initial feasible 
candidate Lyapunov function which is then improved via 
an iterative process. To construct this initial candidate we 
look at the dynamics linearized about the trajectory by 
constructing A : [r, i f ] ^ R nxn , with A(t) = ^(t,x (t)). 
Then, fixing a function Q : [r, if] §™, we solve the 
Lyapunov differential equation: 

-fff(t) - A(tyP*(t)+P*(t)A(t) + Q(t), P *(i f ) = Pf (5) 

over the interval [r, if] . The following Lemma suggests a 
procedure for finding a candidate Lyapunov function. 

Lemma 4- Given a solution to tne Lyapunov equa- 

tion above, there exists a positive constant c such that V* 
defined in equation (4) with: 

P*(t)=exp(c|^)Po*(*) 

satisfies the conditions of Lemma (3), so that: 

T = {(t, x) e [t, if] x M n | V(t, x) < 1} 
is a positive funnel into Q. 

Proof. We begin by_changing coordinates to x(t) = x(t) — 
Xo(t), and defining f(t,x) by: 

x = /(i, x) = f(t, x(t)) - f(t, x (t)) 

We can decompose f(t,x) as: 

x = A(t)x + f(t,x) 

where /(t, x) consists of the second and higher-order terms 
in f(t, x). Taking 

V(t,x) = exp (c^T^) x{t)'P$x{t) 

we have 

V(t, x) = exp (c^— ^ [2x'P^f(t, x) + x'(P* - cP*)x] 

= exp (c^^j [2x'P$f(t, x) - x'(Q + cP*)x}. 

Now, (9fi t/ is a compact set, so f(tf,x) is bounded for 
x e dfl tf , therefore there exists a sufficiently large c = c\ 
such that V(tf,x) < for x e 9fi t/ . 

Since P(t) > and is continuous in i and V is continuous 
in x and i, there exists a time t m <tf such that V(t, x) < 
for all x € (9fit for all t e [i m ,i/]. 

We now show that this is also true on i e [r, i/]. Since 
*m < f° r any e > there exists a c sufficiently large 
that the sets <9fi t are contained in the ball \x\ < e for all 
t E [t, t m ]. 



Since x'P */(i, x) contains third and higher orders in x, 
there exists a sufficiently small e such that \x(t)' Pg(t)f(x, t)\ < 
x(t)'{Q{t) + cP$(t))x(t) for all x <= dn t for all t G [r,i m j. 

This implies that there exists a sufficiently large c = c 2 
such that V"(t, a;) < for all x £ dVL t for all t € [t, t m ]. 

Taking c = max{ci, C2} proves the Lemma. □ 

As we only consider a finite time interval, this result 
guarantees the funnel will have a non-empty intersection 
with {t} x R n for each t € [r, t { ]. 



maX -{ftA,ft}i€Af 

subj. to 



PN-l(t{) 



(8) 



1, 



2.5 Polynomial Lyapunov Functions 



To exploit SOS programming, we develop an alternative 
Lyapunov function to (4) defined by a piecewise polyno- 
mial functions: x : [r, if] i->- M", P : [t, if] h-» §^ and 
p : [r, if] 1 ^ (0,oo): 



V(t,x) = \\x-x Q (t)\\p(t)> p (t) = 



Po(t) 
P(t) ' 



(6) 



In particular, we have Xo(i) approximate x (t) and Po(t) 
approximate Po(t). The function p(t) is a time- varying 
rescaling which we describe in the next section. 

Note that the approximate nature of x (i) docs not pre- 
clude V(t, x) from verifying a funnel as the conditions in 
Lemma 3 are concerned with the level-set where V(t, x) = 
1. As a result, the exact behavior of the function near 
V(t,x) — is not essential. 



3. OPTIMIZATION PROCEDURE 

Once we restrict ourselves to a class of piecewise polyno- 
mial Lyapunov functions, we can approach our optimiza- 
tion task as a bilinear Sum-of-Squares program. In par- 
ticular, we will see that the conditions of Lemma (2) will 
be tests of polynomial negativity on semi-algebraic sets. 
To verify conditions on these specific time-intervals and 
level-sets, we make use of the polynomial 5-procedure (see 
Parrilo (2003)). We arrive in a program with constraints 
bilinear in our parameterization of V(t,x) and multipliers 
used in the 5-procedure. In particular, for this work we 
parameterize V solely by our choice of time- varying p{t). 

Let {ti}fL be a set of knot points which contains the knot 
points of f(t,x), xo(t), and Po(t). In particular, we order 
the knots so that tj < t i+ i for i £ J\f = {0, . . . , N — 1}. 



We define: 



V(t,x) 



x (i)| 



2 

Po(t) 



(7) 



For each interval we define Lagrange multipliers 

£i,Hi € x]_. Let p i} /, and Vi be the polynomial pieces 
of p, f and V on the these intervals. For some constant 
e > 0, we will attempt to optimize a cost function h(p) 
according to: 



VieAf: 

Pi{t) e s[t],pi(t i+ i) = pi + i(t i+ i), 

^Vi(t,x)fi(t,x) + ^Vi(t,x 
- pi(t) + Hi{t,x){pi{t) - Vi(t,x)) 



£ — 



+ £i(t,x)(t-ti)(t i+1 -t) 



e E[t,x] 



£i e Z[t,x}. 

We note that the volume of the set T defined by V(t, x) 
is proportional to: 

vol (J 7 ) cx 




P(t) n 



dct(PoW) 



dt. 



As a surrogate for this cost function, we optimize the linear 
cost: ^ 

h(p) = I p(i) di 



['pit) 



which can be computed exactly. 

We demonstrate that if a given set {pi, (i, Pi}i<£j^ is fea- 
sible, then p(t) defines a function V(t, x) satisfying the 
conditions of Lemma 2. First, we see that p is constrained 
to be continuous and positive so that V(t, x) will be piece- 
wise continuous, piecewise continuously differentiable and 
positive. 



Writing Vi = we have: 



d 



dx 



V i (t,x)f i (t,x) + —V i (t,x) 



dt 



(9) 



equivalent to: 
1 



— V i (t,x)f i (t,x) + -V i {t,x) 



Pi(t) 



Vi(t,x) 
Pi(t) 



As a result, for t e [ti, t i+1 ] and x such that Vi(t, x) = pi{t) 
(equivalently, x e dCl t ) the optimization (8) verifies: 

— Vi{t,x)fi(t,x) + g-Vi{t,x) < £ P {t)-\ 

Note that these constraints verify that: 

^-V(t,x)f(t,x) + ^-V(t,x)<0 

where we impose the constraint on both left and right 
derivatives with respect to time where is not continuous. 

This optimization is unfortunately bilinear in the coeffi- 
cients of p and the Lagrange multipliers. However, this is 
amenable to an alternating search. Given a feasible p(t), 
we can compute Lagrange multipliers. Then, holding these 
multipliers fixed we can attempt to improve p(t). This is 
described in detail below. 



3.1 The Multiplier- Step: Finding LaGrange Multipliers 

For a fixed p(t) we can compute Lagrange multipliers 
via the following set of optimizations. For each interval 
[U,ti + i] we optimize over two polynomials, pi and li in 
R[f, x], and a slack variable jf. 



subj. to 74 



(10) 



— y 4 (t,x)/ l (t,x) + -^(t,x) 



+ ^i(t,a;)(t-ti)(ti + i-t) 



e EM, 



4(i,:r) G E[i,a;]. 

These programs can be computed in parallel. If each 7$ is 
negative, then the set {/Xj, £j, Pi}ieN is a feasible point for 
the optimization (8). 

An obvious question is how to first obtain a feasible p(t). 
Motivated by Section 2, we suggest the following search. 
We search over a positive constant c > 0, and take p(t) to 
be a continuous piecewise polynomial approximation: 



p(t) w cxp -c 



U-t 

tf - T 

If a given choice of c does not verify a funnel (i.e. if the 
optimal values of (10) are not all negative), we iteratively 
increase the value of c, and potentially the number of knot 
point of p(t). 

3.2 The V-Step: Improving p(t). 

Having solved (10) for each time interval, we now attempt 
to increase the size of the region verified by V by opti- 
mizing to increase p(t). To do so, we pose the following 
optimization with {£i, Pi}i^j^ fixed from a solutions to 
(10): 

max -{«b e ^ Hp) ( n ) 
subj. to p N -i{tf) = 1, 
VieAf: 

Pi(t) g E[t],pi(t i+1 ) = p l+1 (t l+1 ), 



e — 



pi(t) + m(t,x){pi{t) - Vi(t,x)) 



+ l i {t,x)(t-t i )(t i+1 -t) 



e Y,[t,x}. 



So long as the class of pi(t) includes the pi(t) used in the 
optimizations (10), this optimization will be feasible, and 
can only improve the achieved value of h(p). 

4. SAMPLING IN TIME 

In terms of computational complexity, the most immediate 
limiting factor to the above approach is performing Sum- 
of-Squares optimization is the degree of t in the polynomi- 
als being constrained. We now discuss an approximation 
to verifying the funnels based on sampling in time. 

The proposed approximation validates the conditions of 
Lemma 3 only at finely sampled points. For each interval 
from the above formulation, we choose a finer 
sampling ti = Tn < r i2 < . . . < r iMi = For each 

such Tij, we define a Lagrange multiplier pij e R[x]. We 
pose the bilinear SOS optimization: 



max -{p^ l3 w Hp) 

subj. to p N -l(tf) = 1, 

Vi G A/",V j G {1,...,MJ : 

Pi(t) g T,[t], pi(t i+1 ) = p i+1 (t i+1 ), 
d - 

-q^ Vi {Tij ,x)fi (Tij , X) 



(12) 



d - 

+ -Q t v i{nj,x)- fain,) 

+ Pij(x)(Pi(l~ij) - Vi(Tij,x)) 



G S[x]. 



This optimization removes all algebraic dependence on t, 
however it does not in general provide an exact certificate. 
One would hope that with sufficiently fine sampling one 
will recover exactness. A partial result to this effect exists. 

Lemma 5. Let Vi : [tj,tj+i] x R n and fi : [£,,f i+ i] x 
R" 1 ^ R™ be continuously diffcrcntiable functions of t, x. 
Further, say that for all t G and x such that 

Vi{t,x) = l,&V(t,x)?0. 

Then, if there exists a t G [t,, tj+i] and x € R" such that 
V(t, x) = 1 and: 

g(t,x) := ^V5(t,a:)/ i (t ) a;) + ^V i (t,a:)><J 

for some positive <5, then there exists an open interval 
around t such that for r in that interval there exists a 
y G K" with: 

^( T >2/) = 1, and 9( T ,y) > 0- 

Proof. As g is continuous in x there exists an r\ > 
such that for all z in B(x,rj) = {z \ \\x — z\\ < 77} 
we have g(t,z) > 5/2. As -j^Vi(t,x) 7^ 0, there exists 
a 2i,Z2 G B(x,rj) such that Vi(t, zi) = 1 + £1 and 
Vi(i, 2^2) = 1 — £2 for some constants £1,62 > 0. As ^Vi 
is continuous, and B(x,rj) is bounded, J^V^ is bounded. 
So there exists an interval around t such that, for any 
r in this interval, Vi{r,z{) > l,T / i(r, 02) < 1. As 5 is 
continuous in i there exists a sub-interval such that for r 
in this interval additionally g(r, z) > for all z G B(x,rj). 
As Vi(t, a;) is continuous, there exists a G (0, 1) such that 
for y = 9zi + (l — 9)z 2 we have V^(r, y) = 1. But, as -B(x, 77) 
is convex, y G B(x, r/) so that g(r, y) > 0. 

Clearly the function Vi we have considered thus far satisfy 
the requirements of this Lemma. The result suggests that 
given a V(t,x) that does not satisfy the conditions of 
Lemma 3, there exists a sufficiently fine uniform sample 
spacing such that V(t, x) will not be in the feasible set of 
(12). This result is partial in that it does not construct 
a sufficient sampling bound from computable quantities, 
and further applies only to a fixed V(t,x), whereas we 
optimize over V with fixed sampling intervals. We are 
currently studying ways of integrating more constructive 
bounds into the optimization (12). 

We use an analogous strategy of bilinear alternation to 
approach (12). The same strategy is used to find an initial 
feasible p(t), and then Lagrange multipliers and p(t) are 
iteratively improved. 



5. NUMERICAL EXPERIMENTS 

We first illustrate the general procedure with a one- 
dimensional polynomial system. Our second example is 
an idealized three degree of freedom satellite model. For 
this second example we compare numerically the proposed 
techniques. 

5.1 A One- Dimensional Example 

We examine a one dimensional time-varying polynomial 
differential equation: 



Optimized Inner Approximation of J 7 * 



-x(t) = f(t, x(t)) = x - \x 2 + 2t - 2 At 



(13) 
[0,1]- 



dt y ~' 'v->--v" - 2 - 
over the interval t £ [—1,1]. Our goal region is Q 
As the system has no control input, we can solve nearly 
exact bounds on the region F* C [—1,1] xi which flows 
into the goal by computing the solutions to (13) with final 
value conditions x(l) — 1 and x(l) = 0. 

To find our inner approximation T C J 7 *, we compute an 
approximate numerical trajectory, xo(t) with final value 
Xq(1) = Xf = 0.5. We take Pf = 4 so that £{ = {x | \x — 
x f\% < 1} = [0,1]- We numerically solve the Lyapunov 
equation (5) to determine P(t). 

We use N = 40 knot points, {tAf =1 , chosen to be the steps 
of a the variable time-step integration of the differential 
equation Lyapunov. We interpolate Xa(t) with a piecewise 
cubic polynomial and P(t) with a piecewise linear function. 
To find our initial candidate Lyapunov function, we begin 
by taking p(t) to be a piecewise linear interpolation of 

exp ^^2~^")> f° r c — 0- Taking c = 4 provides a feasible 
candidate Lyapunov function. This feasible solution is then 
improved by bilinear alternation. Both the initial and 
optimized sets are plotted against the known bounds in 
Figure 1. 

After a single bilinear alternation, a tight region is found. 
Note that the symmetry of the Lyapunov function around 
the trajectory restricts the region being verified. Addi- 
tional trajectories could be used to continue to grow the 
verified region. 

5. 2 Trajectory Stabilization of Satellite Dynamics 

We next evaluate a time-varying positively invariant re- 
gion around a feedback stabilized nominal trajectory. In 
past work, Tedrake et al. (2010), it was demonstrated 
how trajectory optimization and randomized search can 
be combined with such certificates to approximate the 
controllable region for a smooth nonlinear system. 

We examine the stabilization of a nominal trajectory for a 
rigid body floating in space subject to commanded torques. 
The state of the system is given in terms of the angular 
velocity of the body about is principal axes, u> £ K 3 , 
and the Modified Rodriguez parameters, tr e K 3 . Any 
trajectory which excludes full rotations can be represented 
by this projection of the quaternion representation of 
orientation. The kinematic equations are given by: 

er 3 cr 2 ] \ 



)J + 2(70-' - 2 



cr 3 a x 
p 2 er 3 0. 




Fig. 1. The ideal set T and inner approximations calcu- 
lated by the method. Surrounding the nominal tra- 
jectory (solid blue) are time-varying intervals. An 
initial candidate Lyapunov function (red open circle) 
is then improved via the bilinear optimization (solid 
green circle). In this case, a single step of alterna- 
tion provided a certificate tight to the known bounds 
(black stars). Note that the certificate is symmetric 
about the trajectory, and as a result is generally sub- 
optimal. 

The dynamic equations are given by: 

Hu = -(oj x Hoj) + u (15) 

where H = H 1 > is the diagonal, positive definite 
inertia matrix of the system and u £ R 3 is a vector of 
torques. In our example H = diag([5, 3, 2]). The state of 
the system is x = [o j ,lj']' £ R 6 . Together, (14) and (15) 
define controlled dynamics: 

x(t) = f (x(t),u(t)). (16) 

We design a control policy u(t) = Tt(t,x(t)) such that the 
closed loop system, 

±{t) = f{t, X{t)) = f (x(t), 7T(t, X(t))), (17) 

satisfies the assumptions of our method. 

Our goal region is defined by an ellipse centered on the 
origin, described by the positive definite matrix: 



Pa 



'36.1704 



12.1205 







17.4283 



7.2723 







9.8911 



4.8482 



12.1205 



9.1505 







7.2723 



7.3484 







4.8482 



6.2557 



The ellipsoidal region, Q = {x | |ja;||p G < 1} was computed 
numerically as an inner-approximation of the region of 
attraction for the system after feedback stabilization of 
the origin. 

We begin with a nominal command: 



(14) 



«o(*) = -^t(*-5)(< + 5) 



-1 
-1 
1 



Table 1. Runtime comparison for exact and 
sample based approaches. 



Progression of Iterations 





Multiplier Step (sec) 


V Step (sec) 


Sampled 

Exact 


111 112 
5316 5357 


220 220 
1336 1420 



on the interval t £ [0,5]. We compute the solution, xo(t) 
to (16) with u(t) = u (t) and x(5) = 0. Next, we design a 
time-varying LQR controller around the trajectory based 
on the dynamics linearized about the trajectory. Given two 
cost matrices, R £ and Q £ S", we solve the Riccati 
differential equation: 

-S* = A'S* + S*A + Q- S*BR~ 1 B'S*, S*(5) = P f 

where Pf — I.OIPq an d: 



A(t) 



dx 



fo(x (t),u (t)), B(t) 



du 



fo(x (t),u (t)). 



This procedure gives us a time- varying gain matrix: 

K{t) = R^BityS^t) 
The ideal policy is given by: 

TT*{t,x(t))=U (t)-K(t){x-Xo(t)). 

To force Tr(t, x) to be piecewise polynomial in t and poly- 
nomial in x we take a piecewise constant approximation 
of K of K and a piecewise cubic approximation xo(t) of 
x (t). Our control policy is then: 

7r(t, x(t)) = uo(t) - K(t)(x - x (t)). 



%o(t)\\s (t)/p(t) 



We now examine the closed loop dynamics using the 
Lyapunov function: 

V[t,x) = | 

where So(t) is a piecewise linear approximation of S*(t). 
All of the above approximations were piecewise with N = 
49 knot points chosen by a variable time-step integration 
of the Riccati differential equation. 

We compare the time required to compute an exact cer- 
tificate versus the sample based approximation verifying 
necessary conditions at Mi = 10 points in each interval 
[U, ti+i] • In the case of computing exact certificates, the ex- 
pression for ^V(t, x(t)) is degree 7 in t. For both processes, 
we found that the exponentially decaying initial p(t) with 
c = 3 was feasible. Figure 2 shows the progress of two 
iterations using the exact method against the initial p(t). 
The figure plots the volume of O t for each t £ [0,5]. We 
have found that typically very few iterations are required 
for the procedure to converge; indeed, in this case the first 
iteration was already very good. Figure 3 compares the 
results from two iterations of the time-sampled and exact 
methods. In this instance they are nearly identical. Table 
1 compares the run times. The Multiplier Step consisted 
of 48 independent SDPs for the exact method and 480 
independent SDPs for the sample based method. These 
computations were performed on a 12-core 3.33 GHz Intel 
Xeon computer with 24 Gb of RAM. These Multiplier Step 
programs could be trivially parallelized for speedup. 

6. DISCUSSION AND FUTURE WORK 

We now discuss simple extensions and implementation 
details for the method. 
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Fig. 2. The result of successive iterations optimizing the 
volume of verified funnels. After two iterations, the 
methods, the volume ceases to improve. 



Comparison of Optimized Funnels 
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Fig. 3. Comparison of optimized volumes using the exact 
method and time-sampled method. 



6.1 Sampling in Time without Splines 

When applying the time-sampling method described in 
this paper, one is no longer bound to finding a polynomial 
Lyapunov function. As a result, one case use the numeri- 
cally determined xo(t) and Pq (t) instead of interpolations. 
We make use of the splined quantities in both the sampled 
based and exact methods in this work only for the sake of 
comparison. 

6.2 Verifying Families of Funnels 

For several applications, it may prove useful to verify that: 
— V(t,x)f(t,x) + -V(t,x) <0 



for all level-sets in a range V(t, x) € [a, 1], with < a < 1. 
For example: when performing real-time composition of 
preplanned trajectories, it may not be known in advance 
what the goal region is for each trajectory segment. By 
verifying funnels for a whole family of goal regions, the 
best verified funnel can be chosen at execution time. 

6.3 More General Lyapunov Functions 

In this paper work we have restricted ourselves to rescal- 
ings of a time-varying quadratic form, centered near tra- 
jectories. Our first numerical example demonstrated how 
this class can be quite conservative. In principle more 
general polynomial Lyapunov function can be addressed 
with the method presented. Using a richer class of poly- 
nomial Lyapunov functions has proven advantageous for 
time-invariant region-of-attraction analysis (for example, 
see Topcu et al. (2008)). We are currently investigating 
extending the given algorithm for this case. 

6.4 Stability of Limit Cycles 

The authors have also adapted the procedure proposed in 
this paper to verify regions of attraction for limit cycles of 
hybrid systems (see Manchester (2010); Manchester et al. 
(2010)).' 
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